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Background: 


This program of research was initiated with the submission in June 10, 1997 of the 
proposal entitled 

“Wavelet and Multiresolution Analysis for Finite Element Networking Paradigms,” 

The original starting date for the proposal was suggested as July 14, 1997. Nearly 1 1 
months later, on 5/1/98, the project was initiated with an award of $20,000 that was split 
between researchers at the University of Florida and the University of South Carolina. 
The original statement of work in the proposal called for 

“. . . a focused plan of research. . . to extend [the authors] recent contributions 
to multiresolution and wavelet analysis to derive, develop and implement: 

(i) wavelet based methodologies for the compression, transmission, 
decoding, and visualization of three dimensional finite element 
geometry and simulation data in a network environment, 

(ii) methodologies for interactive algorithm monitoring and tracking in 
computational mechanics, and 

(iii) methodologogies for interactive algorithm steering for the 
acceleration of large scale finite element simulations.” 

The original budget for the above statement of work was $252,200. During the 
contractual period, the project monitor Dr. Jerrold Housner moved to a different position 
at NASA Langley Research Center, and the project ceased to be funded. A total of 
$20,000. was awarded in total to the principal investigator at the University of Florida, 
and co-principal investigator at the University of South Carolina towards the research 
project. 

Even with a nearly 1 1 month delay in the disbursement of startup funds, a decrease in 
requested funding of less than 9% of the requested budget, and the removal of the 
contract monitor within NASA Langley, the research team was made unusually strong 
progress towards the completion of the research outlined in this proposal. 

Motivation 

The effective utilization of many large scale numerical simulators requires an interactive 
capability to judge validity of the model chosen as well as the accuracy and efficiency of 
the solution procedures applied. Using their best scientific judgment, the scientist must be 
able to rapidly modify simulations for improvement in the fidelity of simulations, or to 
effect design changes in a collaborative environment. Performing this work on remote 
parallel machines poses special problems which must be resolved for effective utilization 
of these resources in meaningful applications. Historically, a diverse collection of 
parallel computational techniques have been developed for a wide class of 
multiprocessor hardware in order to iteratively solve systems of linear equations 



associated with structural analysis, to approximate solutions of coupled design 
optimization problems, or to obtain time accurate solutions of aerodynamic flows over 
complex aerospace vehicle geometries. However, to enable simultaneous or collective 
use of these tools by different analysts in different laboratories, via virtual environments 
that necessarily operate over networks, it is absolutely imperative that a number of as yet 
unsolved technical barriers must be overcome. 

Summary of Accomplishments: 

The first year of research completed a focused plan of research that 
specifically addressed such technical barriers; 

1 . Originally designed specifically for the Intel Paragon architecture, the tracking and 
steering controller library (Kaulgud, A. and R.C Sharpley. 1995. An Interactive 
Tracking/Steering Library. IM3 Report 95:10. Department of Mathematics, University of 
South Carolina, Columbia, SC (Aug.)) was ported by project personnel to the general 
message passing language MPI (W. Gropp, E. Lusk and A. Skjellmun, USING MPI, 
Portable Parallel Programming with the Message-Passing Interface. The MIT Press, 
Cambridge, 1995.) for use with massively parallel MIMD machines. 

2. Additional improvements were made to network code components of the 
tracking/steering library to permit secure use of the library through firewalls. This was 
tested using a fluid flow in porous media model (for a description, see L. Scott Johnson, 
A. Kaulgud, R.C. Sharpley, R E. Ewing, Z. Leyk, J. Pasciak, M. Celia, and JR. Brannan, 
Integration of Contaminant Transport Simulators on Parallel Machines with a Graphical 
User Interface for Remote Interactive Modeling, in Proceedings of the 1997 Simulation 
Multiconference," Atlanta, April 1997, Soc. for Computer Simulation International, San 
Diego.) over the vBNS network between the University of South Carolina and Texas 
A&M University. 

3. The researchers continued to extend their earlier contributions to multiresolution and 
wavelet analysis to derive, develop and implement wavelet-based methodologies for the 
compression, transmission, decoding and visualization of multidimensional finite 
element logically rectangular geometry and simulation data in a networked environment. 
The wavelet library WV was Anther refined, incorporating encoding algorthims 
including the bitstream encoder of Gao, et al (Z. Gao, A. Andreev and R.C. Sharpley, 
Data Compression and Elementary Encoding of Wavelet Coefficients, EVU Report 97:02, 
Department of Mathematics, University of South Carolina, Columbia, SC (Jan. 97)). 

On-line documentation of WV and a web-based demo of application codes for data 
compression and feature extraction are respectively available at 

http : //www . math . sc . edu/~sj ohnson/wvlib/ 
http : //www . math, sc . edu/~sj ohnson/wvlib/demo/ 



4. The research team extended its methodologies for handling other classes of simulation 
data over network environments. These efforts included the development of 
multiwavelet, divergence free formulations amenable for the compression and 
transmission of flow field simulations as discussed in 

Multiwavelets and Particle Image Velocimetry Methods, A. Kurdila, O. Rediniotis, T. 
Strganac, J. Ko, 40 th Structures, Dynamics, and Materials Conference, AIAA Paper 
Number AIAA-99-1309. 

5. The research team has likewise extended its applications of multiresolution-based 
simulation methods to include reduced order model representations of dynamical systems 
appearing in classical aeroelastic studies. In this paradigm, low dimensional models are 
derived via multiresolution and wavelet methods from classical, high-dimensional FEM 
models of coupled fluid-structural interaction problems. The low order models enable 
the possibility of network simulation and real-time transmission of simulation results to 
remote sites for evaluation by engineers. The technical results are summarized in the 
publication 

Multiresolution Methods for Reduced Order Models for Dynamical Systems, A. Kurdila, 
C. Prazenica, O. Rediniotis, T. Strganac, 40 th Structures, Dynamics, and Materials 
Conference, AIAA Paper Number AIAA-99-1263. 


The details of these latter two contributions discussed in (4) and (5) above are included in 
the technical attachments. 



Appendix (I) 


Multiwavelets and Particle Image Velocimetry Methods, A. Kurdila, O. Rediniotis, T. 
Strganac, J. Ko, 40 th Structures, Dynamics, and Materials Conference, AIAA Paper 
Number AIAA-99-1309. 
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Multiwavelets and Particle Image 
Velocimetry Methods 

Andrew J. Kurdila’ 

Department of Aerospace Engineering, Mechanics, and Engineering Science, 
University of Florida, Gainesville, Florida 32611-6 

Othon Rediniotis*, Thomas Strganac 6 , Jeonghwan Ko f 
Department of Aerospace Engineering, 

Texas A&M University, College Station, Texas 77843-3141 


Abstract 

In recent work the authors multilevel filtering for the 
analysis of digital Particle Image Velocimetry (PIV) 
based on multiresolution analysis. The essential 
contribution of this work was twofold. We sought to 
demonstrate the feasibility and amenability of 
wavelet based methods for local filtering and cross 
correlation calculations required in PIV methods. 
This work focussed on the derivation of windowed 
cross-correlation expressions for wavelet-based 
expansions that are not orthogonal (or biorthogonal) 
over the cross-correlation window. The methodology 
makes use of recently introduced “refinable 
functions” and generalized connection coefficients 
derived in wavelet-based finite element methods. In 
addition, we sought to develop divergence-free bases 
for flow modeling and order reduction. This paper 
extends these recent developments by deriving 
multiwavelet constructions applicable for the analysis 
of PIV. In contrast to our previous work, (1) the 
techniques herein make use of multiple scaling 
functions, (2) the local cross-correlations do not 
require "boundary modifications" that are 
computationally expensive, (3) the derived wavelet- 
based multiresolutions are orthogonal over the 
interrogation window, and (4) the same 
multiwavelets yield constructions of simple 
divergence free bases suitable for order reduction. 

(1) Review and Motivation 

To motivate the derivation of wavelet-based 
PIV algorithms that follow, we briefly review the 
methodology derived in Ref. 26. Frequently, 
instantaneous planar velocity distributions are 

^Associate Professor, Member AIAA 
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l Copyright © 1999 by the authors. Published by the 
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derived from two captured images, separated by a 
known time interval, by cross-correlating 
corresponding sampling windows in the two images. 
In typical approaches 3 , the image data from a window 
(for a 512x512 image, a typical window size is 32x32 
pixels) taken in the first image and that from a 
window at the same position in the second image are 
cross-correlated. Consider the two consecutive 
images. We let f(kj) and g(k,l) represent the pixel 
intensities (typically 0-255) at pixel locations (kj) in 
the first and second images, respectively. From any 
of a number of texts (see for example, Shapiro and 
Rosenfeld 21 ) treating digital image processing, we 
find that the discrete normalized cross-correlation 
R(m,n) associated with discrete functions f(k,l) and 
g(kj) is given by: 

/?(«,«> 


Rijn, ft) 

lETSwiIiTur s/ftf 

Jfc l 

( 1 ) 


In this expression, the function R(m,n) measures the 
correlation of the discrete functions f(kj) and g(kj) 
when they are relatively shifted by (m,n) pixels. We 
do not explicitly express the limits of summation, as 
this will clearly depend on the window size and 
location. The location of the cross-correlation peak 
gives the mean displacement of the particles in the 
interrogation window. This process is repeated by 
moving the window, until the entire image is 
covered. In Ref. 3, inaccuracies are introduced by 
the fact that the window in the second image is at the 
same location as the window in the first image. This 
could yield erroneous predictions since some of the 
particles in the first window could, due to their finite 
velocity, move outside of the interrogation window. 
The likelihood of erroneous predictions of the 
particle displacement grows as the flow velocity or, 
equivalently, the displacement increases. Moreover, 
the successful performance of the cross-correlation 
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requires the use of a rather large interrogation 
window, typically 32x32 pixels. The algorithm 
returns a single velocity prediction for the entire 
interrogation window, therefore, all of the points 
within the window are assumed to have the same 
uniform velocity, and any spatial velocity gradients 
within the window are missed altogether. This 
restriction, can be somewhat alleviated by 
overlapping the interrogation windows. In general, 
for a 5 12x512 image, the velocity map has a grid size 
of 16x16. This resolution is inadequate when fine 
length scales are to be resolved. The determination of 
a local maxima is clearly difficult without significant 
filtering. 

As shown in detail in Ref. 26, it is possible to derive 
an expression for the windowed urt-normalized cross- 
correlation for functions in ID expressed in terms of 
wavelet expansions as 

4(2-^) = zx S 

k l m = 0 

( 2 ) 

where the coefficients are the so-called 

connection coefficients, or refinable integrals. In 
two dimensions, the cross-correlation becomes 


- 0,0 
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(3) 


We note the following : 


(i) The integrals defining T°;° are not trivial to 

calculate in general. For many cases of interest, 
the wavelet scaling functions <j> comprising the 
integrals cannot be expressed in closed form. 

(ii) Numerical methods exist for calculating the 
entries r°;° to any degree of precision. 

Techniques for evaluating these integrals are 
discussed in Ref. 23 and Ref. 24. 

(iii) Once the entries of T^ 0 ; 0 have been calculated 
numerically, they can be stored and applied very 
rapidly. The number of nonzero entries in r°j° is 

proportional to the length of the mask defining 
the wavelet. 

(iv) The entries of r°;°are different for different 
(families of) wavelets. 


The expressions for the wavelet-based local cross- 
correlations in equations (2) and (3) allow local 
wavelet filtering, and the development of multilevel 
methods. As an example, the performance of these 
multilevel algorithms is summarized in Figures (1) - 
(9), which show that the multilevel filtering can be 
effective in improving the fidelity of cross- 
correlation calculations. For example, Figures (1) - 
(4) show that local cross correlations can be 
recursively calculated to improve local maximum 
identification using the multilevel algorithms 
described in Ref. 26. If the simple multilevel 
filtering algorithm (zero-pass coarse, all-pass details) 
depicted in Figure (5) is employed, highly accurate 
velocity fields can be reconstructed as depicted in 
Figures (6)-(9). 

Still, the cross correlation expressions in equations 
(2) and (3) are obviously more complicated to 
implement than the conventional cross-correlation 
depicted in equation (1). Simply put, the complexity 
of the former expressions can be attributed to the fact 
that the selected (Daubechies compactly supported) 
wavelets employed in Ref. 26 are not orthogonal over 
the interrogation window. In this paper, we derive 
new classes of multiwavelets that are indeed 
orthogonal over the interrogation window. With 
these new classes of wavelets and filters, the 
calculations required to implement wavelet-based 
PIV methods are greatly simplified. Tables (1) and 
V2) list the newly derived filters corresponding to the 
orthogonal multiwavelets developed precisely for 
local interrogation windows. 

Finally, it has been shown in Ref. 27 that divergence- 
free wavelet bases derived from the work in Ref. 
(28), (29), (30) can provide the foundation for order 
reduction methods applicable to incompressible 
flows. We conclude this paper by deriving 
divergence-free multiwavelets from the same 
underlying orthonormal multiwavelets used in the 
multilevel filtering operations described above. 

(2) Multiwavelets and Multiresolution Analsysis 

The essential difference between the methodology 
introduced in this paper, and the earlier work by the 
authors in Ref. 26, 27 is that we make use of a 
collection of r real-valued, scalar functions to 
generate the multiresolution. In vector form, the 
scaling functions (or generators) and wavelet 
functions are denoted 
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4 [0*00 1 A IVO) 

0'OOJ [Y r (x) 

Carefully note that the superscript denotes the 
specific generator in a family of generators. The 
translates and dilates of the scaling functions and 
wavelets are given by 

f, t (*)J K»<*>J 

In these equations we have introduced the usual 
shorthand notation 

f it (x) = 2‘ l2 /(2‘x-k) 

Each vector of scaling functions and wavelets 
satisfies a matrix two scale relationship, shown 
below. 

4>(x) = X V2[a. ]d> (2x - k) 

v (.r) = X v/2[^ 4 ]<d (2x - k ) 

k 

s 

5 

In this paper, we will deal with orthonormal families 
of multiwavelets. Specifically, we will employ a 
family of multiwavelets derived in Ref. 35 using the 
intertwining techniques described in Ref. 34. The 
requirements that the multiple scaling functions and 
multiple wavelets defined above are orthonormal 
manifests itself as conditions on the matrix masks. 
The requirement that the generators and wavelets are 
orthonormal to themselves is easily denved from the 
identities 

j <h (x)<f> T (x - k)dx =5 oi I 

ft 

l'¥(x)'¥ r (x-k = I 

ft 

Upon expansion, we have 




f x v2[a (2x - s)x <f> ‘ {2x -2k- myj2[a m f dx 

ft * f" 

= i v2[a ]j«j>( 2x - ,v)<f> 1 (2x -2k- m)dx42[a m ^ 

ms * R 

= - 2 K]5 0 r.k^X a J = 28 M • 1 

WJ 

Hence, the masks defining the generators and the 
wavelets must satisfy 

£ 2 0 2 *-m]K] r = 28 7 

m 

z2[b 2 ,J[bJ = 25 

m 

In addition, we require that the generators and 
wavelets define (orthogonal) complementary spaces. 

That is, the generators are orthogonal to the wavelets 
on any fixed level. 

J 0> (-X)*F r (x- k)dx = 0 

ft 

When we expand these expressions, we obtain 
Jx d2[a s )b{2x- s)x <t> r (2„\*- 2 k~ m)\2[b nt ] ! dx - 0 

R * 

x V2[aJJ®(2.x- j)<I> r (2.v - 2k- m)dxd2\bjf - 0 

xij ft 

Thus, the generators will be orthogonal to the 
wavelets on the same level provided that the masks 
satisfy 

XK..p„r=° 

m 

The specific set of orthonormal multiwavelets 
considered in this paper are depicted in Figures (14) 
and (15). Their matrix masks are summarized in 
Tables (1) and (2). 

For the derivation that follows in the next section, we 
will likewise have occasion to interpret the wavelet 
identities above in terms of the Z transform. If 
[hj n- 2,-l,0,l,2,--- is a sequence of 

matrices, the formal Z transform of the sequence is 
defined to be 

nsZ 

Sometimes we simply refer to this as the frequency 
domain representation of the sequence 
[hj n= 2,-l,0,l,2,---. In any of a number 

of standard texts, it is shown that wavelets can be 
interpreted as a two channel filter bank depicted in 
Figure (16). The two channel filter bank is 
comprised of a cascade of convolution, upsampling 
and downsampling operators. It is not difficult to 
show that the mapping from input to output in the 
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two channel filter bank depicted in Figure (16) can be 
derived as 

y x (z)-k[H{z)X{z)+ H(-z)X(-z)] 
v : (z)= {[G(z)X(z) + Girz)X{-z)] 

x - Hy\ + Gy, 

= \{H{z)H(z)X(z) + H{z)H(-z)X(-z) 
+Gz)G(z)X(z) + G(z)G(-z)X(-z)} 

= + {//(z)tf(z) + G(z)G(z)}X(z) 

+ i{H(z)H(-z)+ G(z)G(-z)}X{-z) 

Conditions that are sufficient to guarantee that the 
input to the filter bank exactly matches the output of 
the filter bank can be seen by inspection to be 

H(z)H(z) + G(z)G(z) = 2/(z~' ) 

H(z)H(-z) + G(z)G(-z) = 0 

These two conditions are referred to as the perfect 
reconstruction, and alias-canceling conditions, 
repsectively. A two channel filter bank that satisfies 
these conditions defines a family of biorthogonal 
multiwavelets and their associated multiresolution 
analyses. 

(3) Simple Divergence-Free Multiwavelets 

As shown in Ref 27, divergence free wavelets can 
be employed in many postprocessing tasks associated 
with particle image velocimetry. During numerical 
experiments carried out in Ref. 27, it was found that 
the treatment of finite domains typically encountered 
in digital PIV image processing proved to be 
particularly troublesome. In this section, we show 
that a very simple modification of the methodology 
suggested in Ref. 33 yields families of divergence 
free multiwavelets. The construction is trivial, in 
fact, if we are given an original orthonormal 
multiwavelet basis system. 

Suppose we are given the orthonormal multiple 
scaling functions and wavelets as discussed in 
Section (2). In this case, it is trivial to show that the 
derivatives of the generators satisfy a two scale 
relationship of their own. 

*_(*)= *'(*)= X 242[a k y{2x - k) 

k 

<& .(*)= X V2[a_ ,]<J> . (2x- s ) 

where the new set of masks are defined to be 


KJ = + 2 [«,] 

Similarly, we can define candidate multi-wavelets as 

(x) - X _ ( 2x - s) 

.> 

where 

[b_J = +2[b,] 

It is important to note that the candidate functions 
0_(x), ¥_(*) above do not constitute an 
orthonormal system of scaling functions and 
generators. In a similar fashion, we can construct 
candidate multiple scaling functions and wavelets by 
integrating the original orthonormal generators and 
wavelets. In this case, the associated two scale 
relations can be written as 

<»,(*)= X v2[a ; ]'J. (2.r- k) 
f.(x)=X V2[* Jo . (2x-k) 

where 

[“..]= =[",] 

[*.,]= Mi 

Again, we emphasize that the multiple set of scaling 
functions (x), ( X ) do not constitute an 

orthonormal set of multi wavelets. However, the two 
sets of functions (0_(x), x F_(x)) and 

(<X> + (x), *F + (x)) do constitute a set of biorthogonal 

multiwavelets. This is easily seen by observing that 
the Z-transform of the masks associated with each of 
these pairs of functions satisfies 

A + (z) = +±A(z) B + (z)= + ±B(z) 

A_(z ) = +2 A(z) B_(z) = +2 B{z) 

As a consequence, we have 

A, G)A \z ') T =(+iA(z))(+2A(-zf) = A(z)Ai-zf = I 

(z-'Y =(+ifl(z))(+2B(-z)' ) = 5(z)B(-r/ = / 
AJz)B.(z- , ) t = (+±A(z))(+2B(-z ) r ) = A(z)B(-z) t = 0 

B. (z)A_(z") T =(+ifl(z))(+2^--) r ) = B(z)A(-z) r - 0 

However, the above set of equations simply 
guarantees that the two channel filter bank 
associated with the differentiated functions 
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(x), 'F, ( x ) and the integrated functions 

<I> + (x), *¥ + (x) satisfy the perfect reconstruction and 
anti-aliasing conditions derived in Section (2). 


specific elements of the r-vectors of generators and 
multiwavelets. 

Conclusions: 


From these observations, it is relatively simple, 
although algebraically tedious, to construct simple 
divergence free bases. While the details exceed the 
scope of this short note, a multiresolution analysis 
comprised of divergence free multiwavelets can be 
constructed in terms of the primal bases 


j-v .(■*>001 

tv (*>.00 j 

j 4> + OM y) } 


¥ 


0,0.2 


-vAx)v(y ) 1 


and their respective dual bases 
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[ 0 ] 
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The reader should carefully note that we have 
suppressed the superscript denoting the specific 
generator (or wavelet) in the list of multiple 
generators (or wavelets). That is, each basis 

appearing on the right hand side of these equations 
should have a superscript that runs between l ... r. 


However, with this notational convention, we can 
succinctly state the multiresolution decomposition of 
divergence free subspaces of the square integrable 
functions as 

Each / e (L 2 (R)f (]{}'■ div/ = 0} 


f ~ i si I 

jnZee E’ fV^AtZxZ 


)v y)} 

(9 ( x *y) 


Again, for notational simplicity, it is understood that 
the above summations are carried out over all the 


In this paper, we have reviewed recent work by the 
authors that derive (I) multilevel filtering methods 
and (ii) reduced order divergence-free bases in terms 
of orthonormal, and biorthogonal wavelets, 
respectively. It was noted that the first set of 
algorithms for multilevel filtering could be expensive 
to implement because the selected wavelet bases are 
not orthogonal over the interrogation window. In the 
second set of algorithms, a biorthogonal basis system 
was employed that is unrelated to the orthonormal 
system employed for multilevel filtering. In this 
paper, we summarize a set of orthonormal 
multiwavelets that have been derived via intertwining 
techniques to address these deficiencies. In short, the 
derived multiwavelet bases (i) are orthogonal over 
the PIV interrogation window, (ii) yield fast, local 
cross-correlation calculations that are amenable to 
multilevel implementations, and (iii) can be used to 
derive simple divergence free bases for 
incompressible flow modeling order reduction. 
Algorithmic implementation and performance are 
discussed in a forthcoming paper. 
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Figure (1) Cross correlation, for multilevel 
algorithms : Denoising on multiple decompositions, 
2 level details, Ref. 26 



Figure (2) Cross correlation, for multilevel 
algorithms : Denoising on multiple decompositions, 
3 level details, Ref. 26 



Figure (3) Cross correlation, for multilevel 
algorithms : Denoising on multiple decompositions, 
4 level details, Ref 26 



Figure (4) Cross correlation, for multilevel 
algorithms : Denoising on multiple decompositions, 
5 level details, Ref. 26 
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Figure (5) Multilevel Decomposition 
and Filtering Schematic, Ref. 26 



Figure (6) Consecutively collected frames, 
velocity reconstruction via multilevel wavelet 
algorithms, t=0 sec. Ref. 26 



Figure (7) Consecutively collected frames, 
velocity reconstruction via multilevel wavelet 
algorithms, t=l sec, Ref 26 



Figure (8) Consecutively collected frames, 
velocity reconstruction via multilevel wavelet 
algorithms, t=2 sec. Ref 26 



Figure (9) Consecutively collected frames, 
velocity reconstruction via multilevel wavelet 
algorithms, t=3 sec, Ref 26 



(a) Original POD Mode 2, 61 x 27 
Figure (10) Proper orthogonal decomposition 
Mode 2, Ref 27 


American Institute of Aeronautics and Astronautics 





(b) First Compression, POD Mode 2, 480 


Figure (11) Proper orthogonal decomposition, 
Mode 2, omit 1 st level details. Ref. 27 



(c) Second Compression, POD Mode 2, 1 53 
Figure (12) Proper orthogonal decomposition, 
Mode 2, omit 1 st + 2 nd level details, Ref. 27 



(d) Third Compression,, POD Mode 2 60 


Figure (13) Proper orthogonal decomposition, 
Mode 2, multiscale filter, omit 1 st +2 nd + 3 rd level 
Ref. 27 



Figure (14) Orthonormal Multiwavelet scaling 
functions, Ref. 35 



Figure (15) Orthonormal Multi wavelet scaling 


functions, Ref. 35 



Figure (16) Prototypical two channel filter bank 
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Abstract 

This paper derives reduced order input-output models 
for a class of nonlinear systems by utilizing wavelet 
approximations of kernels appearing in Volterra series 
representations. While Volterra series representations 
of nonlinear system input/output have been understood 
from a theoretical standpoint for some time, their 
practical use has been limited owing to the 
dimensionality of approximations of the higher order 
(nonlinear) terms. In general, wavelet and 

multiresolution analysis have shown considerable 
promise for the compression of signals, images and, 
most importantly for this paper, some integral 
operators. Unfortunately, causal Volterra series 
representations are expressed in terms of integrals that 
are restricted to products of half-spaces, and there is a 
significant difficulty in deriving wavelets that are 
appropriate for Volterra kernel representations (i.e., that 
are restricted to semi-infinite domains). In addition, it 
is necessary to derive Volterra kernel expansions that 
are consistent with the method of sampling used to 
obtain input and output data. This paper derives 
discrete approximations for truncated Volterra series 
representations in terms of a specific class of 
biorthogonal wavelets. When employing a zero order 
hold for both the input and output signals, it is shown 
that a consistent approximation of the input/output 
system is achieved for a specific choice of biorthogonal 
wavelet families. This family is characterized by the 
fact that all of the wavelets are biorthogonal with 
respect to the characteristic function of the dyadic 
intervals employed to defme the zero order hold. It is 


also simple to show that an arbitrary choice of wavelet 
systems will not, in general, provide a consistent 
approximation for arbitrary input/output mappings. 
Numerical studies of the derived methodologies are 
carried out using experimental pitch-plunge response 
data from the TAMU Nonlinear Aeroelastic Testbed. 

Introduction 

A large collection of methods have been investigated 
for obtaining reduced order representations of linear 
and nonlinear dynamical systems in structural 
mechanics, fluid mechanics, aeroelasticity and control 
theory. These methods include such diverse strategies 
as modal synthesis, Ritz vector reduction, rational 
approximation, Hankel approximation and proper 
orthogonal decomposition (POD). In fluid mechanics, 
the study of the underlying qualitative dynamics of 
various classes of flows using POD has been carried out 
in 24 ’ 26 . In these studies, the essential goal is often to 
study the topological dynamics underlying the more 
complex model 26, 251 24 . Sometimes, the ultimate goal 
is the development of control experiments or 
methodologies . In structural mechanics, Ritz basis 
reduction methods denoted “component mode 
synthesis” have emerged as a distinct discipline [see 
Craig, 1981]. Order reduction methods designed to 
preserve the fidelity of specific measurements, 
observations or performance functionals have been 
studied in linear control theory [see Skelton, 1983]. 
Similarly, some researchers have employed reduced 
basis methods in the simulation and control of the 
nonlinear Navier-Stokes equations 30,27,28 More 
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recently, component mode synthesis methods have been 
studied for a class of open loop simulations of 
aeroelastic systems 40 . 

One common feature of most of these methods is that 
they are not directly amenable to online updating or 
adaptive “response subspace” selection strategies. For 
example, the authors have shown that geometrically 
nonlinear control methods can be extremely effective in 
some applications to nonlinear aeroelastic control 
These strategies rely on an accurate reduced order 
model of the nonlinear open loop dynamics. At the 
same time, the authors have shown that some POD 
reduced basis methods can be extremely effective in 
generating low-dimensional approximations of 
uncontrolled flow. As an example, the authors have 
studied the efficiency of POD methods for obtaining 
reduced order models of synthetic jets in 3S . The 
performance of this approach for studying a particular 
response regime is discussed in detail in 38 . By 
including the four POD modes, three of which are 
depicted in Figures (1) through (3), 99% of the total 
energy of the flow is captured by the reduced order 
Navier-Stokes simulation, as shown in Figure (4). 
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Figure 1: POD Mode 1, Rediniotis, Kurdila, ref(38) | 



We emphasize that this study considers an uncontrolled 
response regime. The difficulty, however, is not 
associated with producing a single reduced order model 
that captures the essential dynamics of a particular 
operating or flow regime. Rather, if the overall goal of 


the model order reduction is to enable control of the 
system, the reduced order model must be accurate over 
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Figure 3: POD Mode 3, Rediniotis.Kurdila, ref(38) 
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Figure 4: %Energy retained in reduced POD basis 
simulation, Rediniotis & Kurdila, ief(38) 


a diverse family of response regimes. In other words, 
because the express purpose of control is to alter 
(usually drastically) the system dynamics, any model 
derived from the response history of the open loop, 
uncontrolled system may be a very poor approximation 
for the closed loop system dynamics. This fact is well 
understood, and well documented, in the control theory 
and linear system theory literature [see for example, 
Skelton 1983 and the references therein]. 
Unfortunately, the problem is compounded in many 
applications to aeroelasticity or fluid mechanics in that 
the governing Navier-Stokes equations are inherently 
nonlinear. Some of the richest theory available for 
treating order reduction problems have been derived in 
the context of linear system theory. 

In this paper, we present a methodology that is directed 
precisely towards achieving efficient reduced order 
representations of a class of nonlinear systems that is 
amenable to adaptive and online control methodologies. 
If we acknowledge that the underlying dynamics are 
nonlinear a priori , a reasonable starting point is to 
choose one of the standard nonlinear system 
parameterizations. Possible choices include Fliess 
functional expansions, Chen series or Volterra series. 
Silva in 43 studies the identification of Volterra series 
from impulse response for aeroelastic applications. We 
will consider the identification of similar systems, but 
for general input/output histories. 
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Volterra Operators and Approximation 

In this paper, we will consider those dynamical systems 
expressed in terms of Volterra integral operators, and 
restrict our attention to single-input/single-output 
systems. The output y(t) can be written formally as the 
infinite sum 

y(0=y i (0+ y 2 (0+ U) 


A W 

V„ = 2 J ''l in Vrf^ vl h^)cR,dr\ 

Thus far, we have shown that the zero order hold for 
both input and output induces the discrete Volterra 
Series in equations (7) and (8). We now show that the 
discrete integrals in equation (8) are amenable to 
wavelet induced multilevel approximatioa Suppose 
that the first order kernel h x is approximated as 


where each term y,(t) is the output of the i m order 
Volterra integral operator For i=l or i=2, we have 

(2) 

y 2 {t) = JL t (3) 


(9) 

P 

where the only assumption on the functions cp y p at this 
point is that they are dual to x J m in the sense that 

IVj.pXj.m-S ,.m ( 10 ) 


where the input is u(t). The theoretical foundations of 
the Volterra series in equation (1), sufficient conditions 
for its convergence and the form for higher order terms 
can be found in 42 . In this paper, we are concerned only 
with the first and second (i=l,i=2) order terms of the 
Volterra series. 

Now, we introduce the characteristic function %(t) of 
the unit interval [0,1] 


The first order term in equation (7) becomes 

.Vu., = I J7 (l KjJVj.p (*))x 

(ii) 

JVa. = II 

* p 

which can be rewritten as 


|1 / e [0,1] 

[0 otherwise 


and its scaled and dilated translates 


(4) 


XjAO=2 J,2 X(2 J t-k) 

_j 2 >n 2 J t-ke [0,1] k2- J <t<(k+ 1)2 •' < 5 ) 

[0 otherwise 

We can obtain a zero-order hold approximation of the 
input by writing 


yU, = I ^l,j.n-k-l U J,t 

y u .„ = 2Kj.,-i u j.»-, 

j 


(12) 


Likewise, if we approximate the second order kernel as 


(^.n) 


(13) 

= IIA J .^ j) <Py,(^>Py,(Tl) 

and substitute this kernel in equation (8), we recover, 
for the second order term in equation (7) 


«(0=I“y.»Xy, t (0 ( 6 ) 

After some tedious algebra, it is shown in 39 that the 
zero order hold approximation of both the first and 
second order kernels of the Volterra series is simply 


y», = ir J^II ^,„<Py,(^>Py,(Tl)) 

y 2J.fl t^ r J [r,j).(n-k-l,n-m-l)^ J.k^ J.m 


yj.n~ ^ Kj.s U j,n-t-\ + ^ ^ (^) 

j=n-l r=/i-l j=n-l 

where 


y 2J,n j.k^ j,m 

yij.» = 'Lh 2 ,Jts.p) U j.»-i-l U J.n-p-l 


( 14 ) 
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Indeed, equations (12) and (14) illustrate that, so long 
as the approximating family is biorthogonal with 
respect to the characteristic functions % jm , we obtain 
the discrete input and output Volterra Series derived in 
equations (7) and (8) by a zero order hold of input and 
output In the next section, we discuss the specific 
selection of a biorthogonal wavelet basis to achieve the 
desired approximation order and order reduction. 


so that 

W j =SpanUf Jl } (21) 

kaZ 

where 

v j, k (*) = 2Jll v ( 2J x - k ) ( 22 > 


Biorthogonal Wavelet Approximations 

In the last section, we showed that a consistent 
approximation of the discrete Volterra input - output 
mapping is achieved when kernel approximates are 
selected to be biorthogonal with respect to 
characteristic functions of dyadic intervals. In this 
section we discuss a class of biorthogonal wavelets that 
satisfy this condition and induce a multilevel 
approximation of the Volterra kernels. Recall that a 
multiresolution analysis is a nested sequence of spaces 

-V^cV,cV,cV t - (15) 


A similar definition holds for the dual wavelet y As 
discussed in 31 , it is possible to define dual wavelets 
and \jr associated with the biorthogonal 

multiresolution analyses [V k } ke2 and that satisfy 

= (23) 

In particular, this can hold if we find masks 
KMM-R} a™ 1 {£} ^ satisfy the two scale 
relationships: 

(p(/) = -/2X<Vp(2f-fc) 


where V 0 is the span of the translates of a fixed function 
<P- 


V 0 = SpanUp(x - k)} (16) 

In this equation, Z is the collection of (signed) integers. 
The remaining spaces in the sequence of equation (15) 
are defined by dilation. We define 

q ? Jt (x) = 2 J,1 (p(2 J x-k ) (17) 

and subsequently 

Vj - S/>art{q>, t } (18) 

keZ 

If we have a second multiresolution \V k } generated 

l J kaZ 

by the function <p, we say that the pair {V t } ti2 and 
jpjj f orm a biorthogonal multiresolution provided 
that 


where 


<p(f) = V2Io/p(2f-fc) 


V(0 = V2Ii*<p(2f-*) 

v(/) = V2I6 t cp(2r-fc) 


(24) 


**=(- 

(25) 

If equations (24) and (25) are satisfied, fast multilevel 
algorithms can be derived that project a function onto 
various nested spaces. If / then / has a 
representation 


/ = (26) 


via equation (18), or as 


(<P jt#*) = J«P^( JC )f^(*) <4c = 5 *» (19) 


/ = !«,. 1 ,tP;-u(0+2:P,-uV J -u(^) (27) 


A wavelet y is a function whose dilates and translates 
span the complement spaces W } defined via 

V h =v j + W j (20) 


from equation (20). Given the coefficients <x jk in 
equation (26), it is easy to show that we can find the 
coefficients a^. u and p >u in equation (27) 
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n 

P >u =XP„- 21 a,.„ 

n 


(28) 


As an example, we obtain 




(29) 




Similarly, if we define the tensor product scaling 
functions 


< I > ‘(^fi) = 9‘(^) ( p'(Tl) 
<t*(^,Tl) = <P(4XP(T1) 

and the tensor product wavelets 

^ l (^fi) = 9(^)V( T l) 
^(^fi) =9^)901) 
'P J (^Ti) = V(^)V(Tl) 

^ r (^Ti) = tp'(^)¥'(Tl) 
^ J '(^fi) = V’(^>P’(fi) 
'P 3 ‘(^fi) = v'(^'(T0 


Equation (28) is the decomposition formula associated 
with the set of biorthogonal wavelets (<p,\y)and ((p,\j/). 
A reconstruction formula can likewise be derived that 
gives the fine scale coefficients in terms of the coarse 
scale. 

a J = X ( a - 2 * a J-U + b m-U P J- U ) ( 3 °) 

These expressions can be used to obtain a multilevel 
representation of the Voltena kernels. Assume that the 
first order Volterra kernel is expressed on the finest 
scale as 

^ ) = X \ , ,9 ,., ($ ) = X « ,.,9 (S ) (31) 

P ' P 

where 

(*+1)2-' 

*,,, = JWX,,6)4 = \ j, (32) 

0 k2-> 

That is, we have identified X JtP =< P JtP - With a 

recursive apphcation of equation (20), a multilevel 
expansion of the first order kernel is obtained. 

A 1 (^) = Xa v ,9,,(^) + XP; 0 ,V,.,(^) 
+XP,.,V,.,£) 

+sp 

: (33) 

+XP,. 

P 

W = Xa*/P A ,G>+2 P„V,,(5) 

p 0 (-J. J 


we obtain the following single scale representation of 
the second order kernel 

# 2 (^-n) = Xa J ,, J) 0. (fJ) (^,Ti) (36) 

The corresponding two-scale expansion of the kernel is 

# 2 (^Ti) = Xa j - u -.. ) ^-,. ( »..,(^n) 

ffl.n 

+xPU<-.-,n.,^(i 1 i) < 37 ) 

m,n 

+xp 3 _ u „, ) n,,»,»(^ T i) 

m.n 

+XP 3 -, ( „.„ ) n,(-..)(^ T D 

m,n 

The development of the multilevel representation of the 
(nonlinear) second order kernel bears close resemblance 
to equation (33). We refer the interested reader to the 
details of the implementation of the first and second 
order kernels in 39 . 

Numerical and Experimental Results 

In a series of papers, the authors have derived, 
implemented and tested reduced order models for the 
prototypical nonlinear aeroelastic system depicted in 
Figure (5) on the next page. This system is comprised 
of a NACA0012 airfoil that is capable of undergoing 
large amplitude response in either the pitch or plunge 
degrees of freedom. Any reasonable representation of 
the response of this system is inherently geometrically 
nonlinear owing to the unique design of the carriage on 
which the airfoil is mounted. The details of the design, 
shown schematically in Figure (6), are described in 
detail in 41 . 
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This system has been designed to exhibit large 
amplitude limit cycle oscillations in particular flow 
regimes, and the performance of geometrically 
nonlinear closed loop control methods is discussed in 
41 . We note that a conventional representation of the 
nonlinear dynamics of this system, that was sufficiently 
accurate to derive closed loop flutter suppression 
control laws, was employed in these previous studies. It 
was comprised of two parametrically dependent, 
coupled nonlinear ordinary differential equations of 
second order. This is, perhaps, the simplest nonlinear 
model that could be used to represent the system. The 
pitch amplitude is so large that it may also be argued 
that stall effects likewise contribute significant 
nonlinear response that is not accounted for in the open 
loop model employed in 41 . In any event, in the current 
numerical example, we study the performance of the 
wavelet-based Volterra series system identification for 
this example. Angle encoders provide a measurement 
of the flap deflection as a function of time, which we 
take as input to the wavelet-based kernel identification 


algorithm. We choose the output to be the pitch angle, 
also measured by angle encoders as depicted in Figure 
(5), measured in radians. In this simple numerical 
experiment, we consider the experimental records from 
a single test in the identification process. The input 
measured flap deflection, output measured pitch and 
flow velocity in the tunnel are depicted in Figures (7), 
(8) and (9), respectively. 
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Figure 7: Experimental input - flap deflection 
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Figure 8: Experimental output - pitch angle 
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Figure 9: Flow velocity in the test section 
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It is important to note that this experimental data 
departs significantly from the framework presented in 
this paper. Specifically, the system is parametrically 
dependent on the flow velocity shown in Figure (9), 
which is clearly not constant. If we note that the flow 
velocity in the figure ranges between 23 and 24 m/s, it 
might be argued that the system is "nearly stationary" 
over the experiment. The implication is, of course, that 
the simple form of the Volterra kernels presented 
earlier are not applicable, strictly speaking. It is 
possible, however, to derive Volterra series expansions 
that are expressed in terms of time-dependent 
kernels. These expansions are applicable for a class of 
non-stationary systems, including the above-described 
system. We anticipate that the identified Volterra 
kernels should indeed vary parametrically with time for 
the system under consideration. The entire data set is 
comprised of 2048 sample points. Using a sliding 
window of 128 sample points, a weakly nonlinear 
Volterra series representation comprised of only the 1“ 
order and 2 nd order kernel was identified for the system. 
For this specific numerical test, the first order kernel 
was comprised of 4 terms, while the second order 
kernel was comprised of a 4 x 4 array. It should be 
noted that the 4 x 4 array is symmetric (see 42 ), so that 
the cardinality of the Volterra series model is 4 + 
(5*4/2) = 14 terms. As shown in Figure (10), graph of 
the predicted output of the Volterra model is not 
discernible from the graph of the experimental output. 



Finally, the evolution of the kernels that comprise this 
non-stationary system can be appreciated by 
considering Figures (11) through (14). Figures (11) 
and (12) depict the representation of the Volterra 
kernels identified from the first and second sample 
windows of the identification process. Each is 
represented in terms of the same 14 basis functions, 
although it is clear that there is a slight variation 


between the sample windows that can be attributed to 
the non-stationary nature of the system under 
consideration. 
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Figure 12: The 2 nd sample window Volterra kernel 


In comparison, the Volterra kernels identified for the 
13 th aril 14 th sample windows, while quite similar to 
one another, vary significantly from the Volterra 
kernels for the 1 st and 2 nd sample windows. Again, this 
is anticipated owing to the non-stationary character of 
the physical dynamical system. 





Figure 13: The 13 th sample window Volterra kernel 
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Conclusions and Future Work 

This paper has derived a wavelet and multiresolution 
based methodology for obtaining reduced order 
approximations of Volterra series. While Volterra 
series representations provide a succinct 
characterization of nonlinear system response in 
principle, their use has been limited in practice due to 
the large number of terms required to represent the 
higher order, nonlinear terms. We show that a 
consistent approximation of the Volterra input/output 
representation is achieved if two conditions are 
satisfied: 

(1) a zero order hold is used for the input and output 
sequences, and 

(2) a biorthogonal wavelet family is selected such that 
the generator is dual to characteristic functions that 
define the zero order hold. 

The identification of a prototypical nonlinear 
aeroelastic system is studied to evaluate the potential of 
the derived method. The prototypical aeroelastic 
system undergoes large amplitude, limit cycle 
oscillations. The experimental data studied in the 
numerical examples is non-stationary, due to the non- 
negligible variations in wind tunnel velocity. 
Nevertheless, the wavelet identification of the nonlinear 
response was extremely accurate with the reduced order 
wavelet models. For a sample record size of 2048 data 
points and sliding window of 128 data points, the 
nonlinear response character was captured with as few 
as 14 wavelets. It was also shown that the nonlinear 
second order kernel did evolve in time, which is to be 
expected for the non-stationary nonlinear model. 
However, the variation in the nonlinear second order 
kernel was essentially slowly varying. One implication 
of the numerical tests is that the migration of these 
methodologies to on-line identification of Volterra 
kernels should be investigated immediately. 


This paper suggests several subsequent lines of 
research. While the ability of the wavelet 
representations to compress integral operators was 
exploited implicitly, the methodology does not 
currently make use of the explicit multilevel structure 
that is available. For example, it is anticipated that 
multilevel and multigrid methods can be used to 
improve the convergence rate of the identification 
procedure. Essentially, the method would perform 
multiscale filtering of the input and output sequences a 
priori, before the kernels of the nonlinear kernel are 
estimated. The kernels themselves could then be 
estimated on resolutions that correspond to the filtered 
input and output sequences. Additionally, this paper 
also suggests that there should be a careful study of 
energy exchange in time-scale space for classical 
nonlinear dynamical systems of aeroelasticity. The 
structure of the kernels may be indicative of the time- 
frequency evolution of energy in the system. Finally, 
the identification procedure suggested in this paper is 
most useful when we can derive associated 
compensators for the nonlinear Volterra series. Such 
work has been studied in some classical texts, but the 
specific forms of the desired compensators for the 
multiresolution kernels have not been studied. 
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